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Abstract 

A new concept is introduced for the adaptive finite element discretization 
of partial differential equations that have a sparsely representable solu- 
tion. Motivated by recent work on compressed sensing, a recursive mesh 
refinement procedure is presented that uses linear programming to find 
a good approximation to the sparse solution on a given refinement level. 
Then only those parts of the mesh are refined that belong to large expan- 
sion coefficients. Error estimates for this procedure are refined and the 
behavior of the procedure is demonstrated via some simple elliptic model 
problems. 
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1 Introduction 

The sparse representation of functions via a linear combination of a small num- 
ber of basic functions has recently received a lot of attention in several mathe- 
matical fields such as approximation theory [18, 33, 39, 38] as well as signal and 
image processing [6, 8, 9, 10, 11, 12, 14, 21, 22, 23, 24, 25, 26, 27]. In terms of 
representations of functions, we can describe the problem as follows. Consider a 
linearly dependent set of n functions 0^, z = 1, 2, . . . , n, (a dictionary [15]) and 
a function / represented as 



Since the set of functions is not linearly independent, this representation is not 
unique and we may want to determine the sparsest representation, i.e., a rep- 
resentation with a maximal number of vanishing coefficients among xi, . . . ,Xn- 
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In the setting of numerical linear algebra, this problem can be formulated as 
follows. Consider a linear system 

= b, (1) 

with 3> G R™'", where m < n and h G W\ The columns of the matrix $ and 
the right hand side b represent the functions and the function /, respectively, 
with respect to some basis of the relevant function space. The problem is then 
to find the sparsest possible solution has as many zero components 

as possible. This optimization problem is in general NP-hard [28, 34]. Starting 
from the work of [14], however, a still growing number of articles have developed 
sufficient conditions that guarantee that an (approximate) sparse solution x 
to (1) can be obtained by solving the linear program 

min ||a;||i, s.t. ^x = b {\\^x - b\\ < e), 

which can be done in polynomial time [30, 31, 32]. We will give a brief survey 
of this theory in Section 2.2. 

In the literature, the development has mostly focused on the construction of 
appropriate coding matrices $ that allow for the sparse representation of a large 
class of functions (signals or images). Furthermore, properties of the columns of 
the matrix (or the dictionary) have been investigated, which guarantee that the 
computation of the sparse solution can be done efficiently via a linear program- 
ming approach, see, for instance, [12, 30]. Often the term compressed sensing 
is used for this approach. 

In this paper we consider a related but different problem. We are interested 
in the numerical solution of partial differential equations 

Lu = /, 

with a differential operator L, to be solved in a domain O c ffi*^ with smooth 
boundary T and appropriate bomidary conditions given on F. 

Considering a classical Galerkin or Petrov-Galerkin finite element approach, 
see e.g. [5], one seeks a solution u in some function space U (which is spanned 
by 01, ... , (pn), represented as 

n 

u = ^Ui4>i. (2) 
1=1 

Again we are interested in sparse representations with a maximal number of 
vanishing coefficients ui. In contrast to the cases discussed before, here we 
would like to construct the space U and the basis functions (jii in the finite 
element discretization in such a way that first of all a sparse representation 
of the solution to (2) exists and second that it can be determined efficiently. 
Furthermore, it would be ideal if the functions 0i could be constructed in a 
multilevel or adaptive way. 

The usual approach to achieve this goal is to use local a posteriori error 
estimation to determine where a refinement, i.e., the addition of further basis 
functions is necessary. For example, in the dual weighted residual approach [3] 
this is done by solving an optimization problem for the error. 
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In this article, wc examine the possibiHty to use similar approaches as those 
used in compressed sensing, i.e., to use £i-minimization and linear program- 
ming to perform the adaptive refinement in the finite element method in such 
a way that the solution is sparsely represented by a linear combination of basis 
functions. In order to achieve this goal, we propose the following framework. 

We determine u e U as the solution of the weak formulation 



Here, V is a space of test functions and (•,•) is an appropriate inner product. 
In the simplest version of a two-level approach, we construct finite dimensional 
spaces of coarse and fine basis functions U" C C U and corresponding 
spaces for coarse and fine test functions V" C C V. Then we determine the 
sparsest solution in , such that 



via the solution of an underdetermined system of the form (1). Based on the 
sparse solution, we determine new coarse and fine spaces U2 C C U, V2 C 
Y2 C V, and iterate this procedure. 

This framework combines the ideas developed in compressed sensing with 
well-known concepts arising in adaptive and multilevel finite element methods 
[17]. But instead of using local and global error estimates to obtain error in- 
dicators by which the grid refinement is controlled, here the solution of the 
^i-minimization problem is used to c;ontrol the grid refinement and adaptivity. 

Many issues of this approach have, however, not yet been resolved, in par- 
ticular, the theoretical analysis of this approach (see Section 4). We see the 
following potential advantages and disadvantages of this framework. On the 
positive side, the ^i-minimization approach allows for an easy automation. We 
will demonstrate this with some numerical examples in Section 5. On the down- 
side, the analysis of the approach seems to be hard even for classical elliptic 
problems, see Section 4 and due to the potentially high complexity of the linear 
programming methods this approach will only be successful, if the procedure 
needs only a few levels and a small sparse representation of the solution exists, 
see Section 5. 

2 Notation and Preliminaries 
2.1 Notation 

For TO, n € N = {1, 2, . . . }, we denote by R™'" the set of real m x n matrices, 
and by /„ the n x n identity matrix. Furthermore, we denote the Euclidean 
inner product on K" by (•, •), i.e., for x,y G K", 



{v, Lu- f) = for all v eW. 



{v, Lu- f)=0 for all vgW]^\Y, 



n 



1 



n 




For 1 < p < 00, the i!p-norm of x G K' 



n 



is defined by 
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with the special case 

||x||oo := . max \xj\, 

je{l,...,n} 

if p = oo . 

The definition of ||-||p can also be formally extended to the case that < 
p < 1. For < p < 1, ||-||p is only a quasi- norm, since the triangle inequality 
is not satisfied, but still a generalized triangle inequality holds, i.e., for every 

x,y gR" one has 

lk + j/ll^< INI^ + lbll^- 

Finally, for p = and a; € M", we introduce the notation 

||.x||o := #supp(a;), 

where supp(a;) := {j £ {1, . . . ,n} : Xj 0} is the support of x. Hence, ||a;||o 
counts the number of nonzero entries of x. Note that in this case even the 

homogeneity is violated, since for a ^ we have HaxHo — \\x\\a. 

For a symmetric positive definite matrix A = G R"'", we introduce the 
energy inner product 

{u, v)a ■■= {u, Av) 

and the induced energy norm 

\\x\\a ■■= \/{x,x)a- 

Every symmetric positive definite matrix A € M"'" has a unique symmetric 
positive definite square root B := As, with A = B"^ = B'^B satisfying the 
relation [29]: 

\\x\\A = \\Bxh. 

2.2 Sparse Representation and Compressed Sensing 

In this part wc survey some recent rcsiilts on sparse representations of functions 
based on the solution of underdetermined linear systems via £i-minimization. 
We also discuss the recently introduced concept of compressed sensing. 

Definition 2.1 ([10]). Let $ e M™'" with m < n and k € {1, . . . ,n}. The 
fc-restricted isometry constant is the smallest number 5k, such that 

{l-5k)\\x\\l<\\^x\\l<{l + 5u)\\x\\l (3) 

for all X E R" with ||x||o < k. 

If $ in Definition 2.1 is orthonormal, then clearly 6k = for all k. Conversely, 
if the constant i5fe is close to for some matrix every set of columns of $ of 
cardinality less than or equal to k behaves like an orthonormal system. In the 
case that < 6k < V^— 1 for large enough k, we say that the matrix $ has the 
restricted isometry property [7, 20] . 

For $ G M"*'" with m < n, a vector of the form b = $.x represents (encodes) 
the vector x in terms of the columns of To extract the information about x 
that b contains, we use a decoder A : R™ — > M" which is a (not necessarily 
linear) mapping. Then y := A(6) = A{^x) is our approximation to x from the 
information given in b. In general, for a given b and matrix $, A(6) may not be 
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unique and it could be a set of vectors. But here for simplicity we take one of 
them and deal with this vector only. 

Let Sfc := {2 G R" : ||z||o < k} denote the vectors in M" of support less 
than or equal to fc. In the following wo use the classical i'p-norm, but also other 
norms are possible, see Theorem 2.2 below. We introduce the distance 

(^k{x)p ■■= min||a;- 2:||j„ 

and observe that for x, 2; e M" and p> 1 the following inequality holds: 

(J2k{x + z)p < ak{x)p + ak{z)p. (4) 
We have the following theorem. 

Theorem 2.2 ([18]). Consider a matrix <I> € R"*'" with m < n, a value k € 
{1, . . . , n}, and let Af = ker(<I>). // there exists a constant Co such that 

Mp < ^ <^2k{v)p, for all v€J\f, (5) 
then there exists a decoder A such that 

\\x- A{^x)\\p<Coak{x)p, forallxGW. (6) 
Conversely, if there exists a decoder A such that (6) holds, then 

WvWp < Co <T2k{v)p, for all ijeM. (7) 

If we combine Theorem 2.2 for p = 1 with the restricted isometry prop- 
erty (3), then we have the following. 

Theorem 2.3 ([18]). Let 3> G R™'", m < n and k € {1, . . . ,n}. Assume that $ 
satisfies 

{l-S3k)\\x\\l < \\^x\\l < il + S3k)\\x\\l 
for all x with ||a;||o < ik, such that 

Ozk<0 < . 

Define a decoder A for $ via 

A(6) := argmin{||a;||i : b = $a;}. 

Then 

\\x - A($a;)||i < Co(Tk{x)i, 

where 

^ V2+1-(V2-1)^ 

Theorem 2.3 shows that the £i-norm solution can be as good as the best 
fc-term approximation. An analogous result is the following. 
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Theorem 2.4 ([7]). Let $ e K™'", m <n and k £ {1, . . . , n}. Assume that $ 
satisfies the restricted isornetry property (3) of order 2k such that S^k < — 1 
and b = + e where ||e||2 < e. // 

A(6) =argmin{||z||i : ||6 - $z||2 < e}, 

then 

\\x^A{b)h<C,^^+C,e 
for some constants C\ and C2 only depending on 52k ■ 

Remark 2.5. It is easy to sec that for the case where e = and x is fc-sparse, 
we have the exact recovery, in other words, x = A(6); sec [7] for details. 

Besides the fc-restricted isornetry constant 5k, a second quantity plays an 
important role in compressed sensing [26, 41, 42]. 

Definition 2.6. Let $ e M™'" with m < n have unit norm columns, i.e., 
$ = [4>i ■ ■ ■ 4'n] with \\(t)i\\2 = 1, for i = 1, . . . ,n. Then the mutual incoherence 
of the matrix $ is defined by 

M{^) :=max\{(t)i,(t)j)\. 
in 

The mutual incoherence M{^) of a matrix i> is related to the fc-restricted 
isometry constant via 

4< (fc-l)M($). 

The following Lemma shows how the mutual incoherence may be used to 
bound the norm of the encoded vector 6 = $a;. 

Lemma 2.7. Let $ = [<j)\ - ■ ■ <f)n] € R"*'" with m < n have unit norm columns. 
Then for every a; e M" the inequality 

\\^xf2<{l-M{^))\\xf2 + M{mx\\l 

holds. 

Proof. The proof follows by the following (in)equalities. 

n n 

= ^^XiXj{(j)i,(j)j) = \\x\\l + ^XiXj{(j)i,(j)j) 

i=l j=l i^j 

< \\x\\l + M($) Y^\x,\ \x,\ = \\x\\l + M{^){\\x\\l - \\x\\l). 

□ 

Lemma 2.7 states that is bounded by a convex combination of ||a;||2 

and with the mutual incoherence as a parameter. 

Compressed Sensing (Compressive Sampling) refers to a problem of "effi- 
cient" recovery of an unknown vector x G R" from the partial information 
provided by linear measurements {x,(f)j),(f)j G K^jj = l,...,m. The goal in 
compressed sensing is to design an algorithm that approximates x from the in- 
formation b = {{x, . . . ,{x, (i)m)) G K™. Clearly the most important case is 
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when the number of measurements m is much smaller than n. The criicial step 
for this to work, is to build a set of sensing vectors 0j G R", j = 1, . . . , m that 
is "good" for the approximation of all vectors x e M". Clearly, the terms "efE- 
cient" and "good" should be clarified in a mathematical setting of the problem. 

A natural variant of this setting, and this is the approach that is discussed 
here, uses the concept of sparsity. The problem can then be stated as follows. 
For given integers m, < n we want to determine the largest sparsity k{m^n) 
such that there exists a set of vectors G = 1, . . . ,m and an efficient 

decoder A, mapping h into R" in such a way that for any x of sparsity k{m, n) 
one has exact recovery A (a;) = x, (see [21]). 

3 Sparse Representations of Solutions of PDEs 

As discussed in the introduction, we want to use similar ideas as those used in 
compressed sensing in the context of the solution of partial differential equations. 

3.1 General Setup 

For a Hilbert space of functions El = IHl(f2) on a domain O c M"^ with smooth 
boundary F, we denote by (/, 5) the inner product of /, g e H and by ||/||h := 
^ (/, /) the induced norm. 

A generating system or dictionary for H is a family {<;ii}^j of unit norm 
elements (i.e., ll(/>i|lH = 1) in H, such that finite linear combinations of the 
elements {4n} are dense in H. A smallest possible dictionary is a basis of H, 
while the other dictionaries are redundant families of elements. Elements of H 
do not have unique representations as a linear sums of redundant dictionary 
elements. 

We will consider elliptic boundary value problems, see, e.g., [5], and we want 
to find the solution of 

Lu = f in Q, 

for a differential operator L and homogeneous Dirichlet boundary conditions 
u = on the boundary F of 

To fix notation and explain our approach, we will lay out a finite element 
approach, where we assume that test and solution space U C H are the same. 
In order to get a closer analogy between the linear algebra formulation and the 
function space formulation, we assume that we have a redundant dictionary 
V = such that 

span{^j-}°ii =U. 

The corresponding weak formulation for the PDE problem is to find u gU such 
that 

a{u, v) := {Lu, v) = (/, v), for all v gV, (8) 

where «(•,•) is a bilinear form. 

Finitely expressing u, v in terms of the dictionary as 

00 00 
u = ^^Ui4>i and v = y^^Vi(pi, 

i=l i=l 
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we can write the problem in terms of infinite vectors and matrices as 

oo 

^ Vj ai, Ui = = (9) 

Here, A°° := [a^] is the stiffness matrix, with := {L(f)i,(j)j), i,j G N, the 
right hand side 6°° is defined by 6?° := {f,(t)i), for i G N, and the coordinate 
vectors are 



Ul 



















The weak solution then can be formulated as the (infinite) linear system 

A°°u°° = b°°. (11) 

Note that if the (f)/s are not linearly independent, the (infinite) linear system (9) 
and hence (11) may not be solvable or may have many solutions. 



3.2 Algorithmic Framework 

Let us now consider finite dimensional subproblems of (11) by assuming that the 
dictionary T> subsumes a refinement procedure, i.e., for each basis function 
there exists a set of refined basis functions in V. In particular, wc assume that 
there is a mapping of each basis function (/>j to the index set ref (i) C N of refined 
basis functions. In a hierarchical refinement, every t^j can be written as a linear 
combination of {(/)j}jgref(i)- 

In many practical applications, the refinement will arise from a geometric 
refinement, for instance, by subdividing some triangulation and corresponding 
basis functions. Furthermore, the refinement may satisfy Tcf{j) nref(^) = for 
j ^ i (but see Section 5). For notational convenience, we define 

ref(5) := |J ref(i), 

for S CN. With this notation, we obtain a sequence of index sets 

T° := {1}, = ref(T°), = ref(T^), . . . 

Define S'^ := L) ■ ■ ■ L) T'^ and denote the corresponding nested subspaces as 

U° C C C • • • C U with U'= := span{0j}jes. . 

We will appropriately select subsets R'^ Q C'^ <Z S^, where i?*^ and are 
seen as subsets of the rows and columns of ^4°° , respectively. The corresponding 
submatrix is defined as follows. 

:=^oo[^fc^C''=] := [(i0i,<^,•)]iei^^iec^• 
The corresponding right hand side is 



8 




Figure 1: Illustration of the stepwise refinement method. Le/t; Picture of a tree that is 
iteratively refined. Right: A view on the matrix ^4°° and its partition as in (13). 



We thus arrive at the finite dimensional subsystem 
and the approximate solution in this case is 



(12) 



jec" 



As in classical adaptive methods the hope is to keep the size of the matrix A'^ 
small (i.e., keep B!^ and C*^ small) and still obtain a good approximation of the 
solution arising from the full refinement of level k, i.e, a solution obtained from 
solving (12) for = = . 

We will now explain how compressed sensing can be used to select small R!^ 
and under the condition that we still obtain a convergent method. For 
this, assume that we have already selected R^^^ and C*^"^. We may start with 
ijo = = {1}, but in practice one should choose an appropriately fine level. 
We now refine these sets to F& := ref(i?'^~^) and := ref(C'^~^). Then A'^ 
and can be partitioned as follows 



A" = 
where 

A21 = A°°[^^C'=-l] 
A31 = A°°fR^C'=-l] 



with R 



All 


A12 


Ai3 


A21 


A22 


A23 


A31 


A32 


A33 



6*= 



A12 = A°°[i?'=-SC*=] 

A22 = A°°[^^c''=] 

A32 = A°°[s^c''=] 



(13) 



Ai3 = A°°[i?'=-\C ] 
A23 = A°°[^^C*'] 
A33 = A°°[e',c'], 



N \ (i?*^ ^ U .R*^) and C " defined analogously. Similarly, the right 
hand side is defined as 

61 = 6°°[i?'=-i], b2 = b°°[R% b3 = b'^\R\ 

The main idea in the construction of R'' and C*^ is to use the underdeter- 
mined subsystem 

[A21 A22] 2; = 62, 
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Algorithm 1 Iteratively Refinement Basis Pursuit (IRBP) 

1: Set i?o = (7° = {1} 

2; for k = 1, . . . , until convergence do 

3: Construct R'' = ref(i?'=-i), C"^ = ref(C"=-i) 

4: Construct A21 = C'"-^], A22 = C% and 62 = 

5; Solve the following minimization problem: 

z'' = argmin{||2r||i : [A21 ^22]^ = 62}. 

6: Let C C C'^"^ U C'^ be the index set corresponding to the support of z'^. 
7: Set C"= = C U C"=-i, R'' := C"=. 
8: end for 



and to compute a sparse solution, by taking a minimal ^i-solution, i.e., 

z'' = argmin{||z||i : [A21 A22] z = 62}- 

For the noisy case, we solve the following problem: 

z'' = argmin{||2||i : ||[A2i A22] z - 62II2 < Cfc}- 

where Efe is a given upper bound on the size of the noise. The submatrix 
[A21 A22] is chosen because it combines the refined rows with the full set of 
columns that are available at the current iteration. 

Now assume that C C C*^"^ U C'^ is the index set corresponding to the 
support of z'^ as defined in Section 2.1. Then the new sets are set to 

C'==CuC'=-\ R'' = C''. 

Thus, the support of z'^ is only used to select basis functions among C'^ and the 
information in C*' H C is not used. 
Note that by construction 

C'' C tree(C), 

where trcc(C) is the set of basis functions on the path of a basis function to the 

root of the rcfincmcnt-tree, i.e., 

tree(i) := e N : 3ji, ...,js with ji € ref(£), € ref(ii), j€ ref(js)}. 

The process is terminated if \\z^ — z'^^^j|2 < e, where e is a given tolerance. 
Since we are using £i-minimization, which in [13] is called basis pursuit, we call 
this process Iteratively Refined Basis Pursuit. The method is summarized as 
Algorithm 1. Figure 1 gives an illustration of the process. 

Example 3.1. By definition, the first sets are = C° = {1}. Now suppose 
that ref (1) = {2, 3, 4}, i.e, the initial rows and columns are = R} = {2, 3, 4}. 
We then solve the ^i-minimization problem for the corresponding 3x4 matrix. 
Now suppose that supp(z-^) = {1,2}, then = {1,2}. Assume that ref(2) = 
{5, 6, 7}, see Figure 2 for an illustration. Then C*^ = {3, 4, 5, 6, 7} and the next 
matrix is of size 5x7, since R^ = {1, 2} and R? = {3, 4, 5, 6, 7}. 
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1 




Note that if the support of is full, then the above procedure yields a 
simple refinement process in which = S''~^ and = . 

In general, this process need not converge. In fact, if at some level we see 
that the error does not decrease, then we back up in the tree of refined basis 
functions and refine at higher levels until we obtain a decrease in the error. 

Remark 3.2. Note that although our description was based on the assumption 
that test and solution space are the same, the principle of the process does not 
depend on this assumption. Similar concepts for adaptive refinement methods 
in the context of wavelets are presented, e.g. in [2, 16, 17]. 



3.3 Properties of the Proposed Method 

In our approach we want to achieve several goals. The solution should be 
sparse, and z'^ should be a good approximation of the solution x'"*"^ of (12). In 
order to analyze the behavior of the suggested approach, we study the case of 
two levels and assume that R'^ = C''. We will also slightly change notation as 
follows. For fc e N, wc set n := # R^-^ and iV := # (R''-^ U rcf Wo 
then introduce the following notation for submatrices of A°° and subvectors of 
6°° as in (13): 



All 

A21 



A12 
A22 



(14) 



and A" := An, 6" := 61. Note that A^ is of size N x N and A" of size nxn. 
The corresponding linear systems are 



A°°x°° = b°°, 



(15) 
(16) 



and j4"x" = 6". In Algorithm 1, we take the last N — n rows of (14) and 
consider the linear system 



\N-n,N „N 



Z =02 



with 



\N-n,N 



:= [^21^22]- 



(17) 



We find a solution of this underdetermined system, by taking a minimal £1- 
solution, i.e., by solving 



argmin{||2;||i : A^ 



(18) 
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As a first step of the analysis we estimate tlie energy norm error between x 



N 



and 



x°°\\ . To derive such a bound, we need 



to embed z'^ and a;" into IR°° by appending as follows: 



zN 











(19) 



We assume that A°°x°° = b°°, A^x'^ = 6^, and that z'^ is determined by (18). 

We want to find necessary and sufficient conditions on the matrix A°° and 
thus on the dictionary such that there exists a constant for which the 
following inequality holds, 



{z'' -x'') ' A'^{z'' -x'' )<C^ ( 



(20) 



The constant should only depend on the ratio Considering the results 
of [12, 18, 26], we may expect that if the matrices A°° , A^ have a small mu- 
tual incoherence or some restricted isometry property, such conditions can be 
obtained. We will come back to this point in Section 4. 

If Inequality (20) holds, then by using the triangle inequality we obtain the 
error estimate 

{x°° - z^fA°°{x°° - z^) < (1 + C^)ix°° - x^fA°°{x°° - xn), 

which means that the (hopefully sparse) solution z^ obtained by solving (18) 
is as almost good as the solution of (16). 

To estimate the errors between x°°, x^ , and z^, we need to consider a refined 
partition of A°° and h°° as defined in (14). In order to relate the solutions at 
different levels of refinement and the solution of the f i-minimization, we make 
use of the following Lemmas. 



Lemma 3.3. Let 



and z^ 



he solutions of the problems (15), (16), and (18), respectively. Furthermore, 
let x^ and z^ he as in (19). Then inequality (20) can he rewritten as 



«-^f)T(6i- [^11^12] 
Proof. Since 



)<Cn (a;r)^(&3-[^3i^32] 



A°°x° 



it follows that 



bi 
b2 

[A31 A32]a;^ 



00 i,N 



A^\x^ -x°^) 



A°°z 



'N' 
Xi 

J-2 



N' 



[AiiAi2]z 
b2 

[^31^32]^^ 



[^31 As^jz"" - 63, 
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and 



Thus, we have 







(f^ - X^fA°°{x^ -X°°) = {X^f{b3 - [^31 



Xi 

X2 



and 



(5^ - - = « - - [An A12] 



^1 



Plugging these expressions into (20) yields the claim. 



□ 



Remark 3.4. A weaker version of (20) and of Lemma 3.3 will be given in 
Section 4. 

The following Lemma gives a condition that has to be satisfied in order to 
guarantee that the refinement process can be iterated. 



Lemma 3.5. Let 



0^ = 



be a solution of (18) , where A^ "'^ is as defined in (17), and suppose that A22 
is invertible. If ^ 0, then 



Proof. Since 



||A-2%l||l>l. 





^22^ 62 



(21) 



is a feasible solution of (18), it follows that 



i<ll^'lli- 



(22) 



Moreover, from ^121-21^ + ^22-22^ = ^2, we obtain that Z2 = A22^^2 — A22^^2i-2i^- 
Thus, using (22), we obtain 

m2-2%||l>||^^||l 

= \\z^h + \\z^h 

= ll^22^''2 — ^22^^212^1^111 + H^^l^Hl 

^ 11^22^ ''2II1 — 11^22^ ^21-2ri|l + ll-Zl^lll- 

Since ||A22^A2izf ||i < 11^22^ A21II1 • ||zf ||i, we have 

INfl|l<IIA-2'^2l||l-||2f||l, 

which completes the proof. □ 
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Remark 3.6. Lemma 3.5 implies that a solution of the i'l-minimization prob- 
lem can only lead to an improvement if the matrix [A21 A22] satisfies (21). 
Otherwise, an optimal solution can already be obtained by solving the linear 
system ^22^2^ = ^2. Another observation is that Lemma 3.5 remains true for 
any nonsingular principal submatrix of . 

In order to compare sparse and non-sparse solutions we introduce the short 
notation s{x) := supp(a;) for a vector x G M™ and s{x) = {1,2,..., to} \ s{x). 
For x e and 5 C {1, . . . , m}, we denote 

{ys)i = <^ . i = l,...,m. 

1 otherwise. 

We then have the following Lemma. 

Lemma 3.7. Let be a solution of (18), where A^~"'-^ is as in (17), and let 
x^ he a solution of A^x^ = , with A^ as in (14). Then for the difference 
§N ._ _ fig^yg ffig inequality 

ll^«V)lli 1 
ll^^lli 

Proof. Since x^ is a feasible solution of (18), we have 

+ = < (23) 

Furthermore, x^ + 6^ = x^ + ^^{x") ^^x'^')■ Therefore, by (23), we have 
lk^||i>||a:^ + 5^||i 

>lk"^l|l-||^«V)lll + ll'^^x-)lll- 

Rewriting yields that 

ll<5,V)lli > = 11'^'" - -^.V)!!! > ll'J^lli - ll<5.V)lli, 

which implies the assertion. □ 
Remark 3.8. The proof of Lemma 3.7 shows that 

||(^^-X%(,.)||l<||(^^-X^),(,«)||l. 

In particular, if {z^ — x'^)s(^x'^) = 0) then we conclude that z^ = x^ . If instead 
of ^i-minimization, we use fo-minimization and compute 

= argmin{||^||o : A^""'^^ = 62}, 

then we get the analogous estimate 

||(«;^-x^)^(,«)||o<||(«;^-x^).(,«)||o. 
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Remark 3.9. In general, it is not true that and a;^ satisfy the inequahty 
||-Ziv||o < ||a;jv||o- For example if 



then 



^3,4 



" 1 


1 


1 


2 


1 1 


1 


V2 2 







1 

2 


1 


1 
2 


1 1 


1 







1 
2 


— l]""" and 


= [0 







z4 = [V2,6,0,-l]T, 
and in this case ||-ZAr||i = 7 + ^/2 < 7 + 2 = ||a;jv||i and 

||2^jv||o = 3 > 2 = ||a;jv||o- 

In this section we have set the stage for the solution of PDEs via ^i-mini- 
mization. In the next section we provide details. 

4 The Restricted Isometry Property for Elliptic 
PDEs 

In this section, we again discuss the special case that solution space and test 
space are the same, i.e., we assume U = V C H and we consider the symmetric 
bilinear form a(-, •) : U x V ^ M associated with the operator L as in (8). 
We also assume that there exist constants ai > 0, a2 < oo, such that: 

Q!i ||w||h < < Q!2 ||w||hi (24) 

i.e., a(-,-) is uniformly elliptic with constant a\ and uniformly hounded with 
constant a2- 

In order to connect this classical norm equivalence with the fc-restricted 
isometry property, we assume that for the dictionary V = the following 

k-equivalence between ||-||e and the £2-iiorm ||-||2 holds, i.e., we assume that 
there exist constants /3i > 0, /32 < oo with 



Pi \\u" 



< 



OO 

IE 

i=l 



U^ <P^\L < h ||m°°||2, 



(25) 



for all infinite vectors u°° as in (10) with the property that ||w°°||o < k. 
Note that inequality (24) can be written as 



< |M°°2 yOO||2 



i=l 



2 < "2 ||E"*'^*IIh' 
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or equivalently as 



(26) 



i=l 



i=l 



Combining inequalities (25) and (26), we obtain 

for all vectors with ||u°°||o < k. 

We consider the dictionary V — {(l)i,(f>2, ■ ■ ■}, and we choose a set = 
{gi, . . . , qn} C N with associated elements (bq^ , • • • , (Ag^ ^ F'^'' tli(3 theoretical 
analysis we may assimie w.l.o.g. that X"^ = {1. 2, . . . , N}. This selection can be 
obtained via an appropriate reordering of the elements 4>i of the dictionary. 
The corresponding finite stiffness matrix associated with this subset is then 



with fly = a{(p,,(l)j), i,j = 1, 



, N. Since we have assumed 



uniform ellipticity and since test and solution space are equal, it follows that 
is symmetric and positive semidefinite with rank(A''^) = iV; is positive 
definite if (f)i, . . . , (pN form a basis. 

Since is symmetric positive semidefinite, A^ has a factorization A^ = 
{B'^)^ , where e K."'^ has full row rank. Hence, there exists a permu- 
tation matrix P such that 



^N-n 



with e M"'" invertible. This yields 
with stiffness matrix A^ = An. Then we have 



An Ai2 
A21 A22 



{PB'^yPB'^ = 





T 


P" 


^N-n 




^N-n 



Suppose that it is possible to choose the permutation matrix P in such a way 
that (measured in spectral norm) 



IP 



N-n\\ 



< e 



with a small e > 0, i.e.. 



||P^-"a;||i<e||a;||i 
for all a; G R". Suppose further that 



(1 - Sk) Wx^'wl < (x^)^^ V < (1 + 4) Wx^'Wl 



or equivalently 



{i-Sk)\\x^r2<\\B^x^r2<{i+Sk) wx^u 



for all x^ with ||a;''^||o < k. 

To get an error estimate between the solution that is based on ^i-minimiz- 
ation and the best fc-term approximation, we first prove the following result. 
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Theorem 4.1. Let A E M.^'^ be symmetrie positive semi-definite, and consider 
the solvable linear system Ax = b. Let A = B^B be a full rank factorization, 
and let P G R'^'^ be a permutation matrix such that the following properties 
hold. 



1. PB = 



B2 



A = B'^B = BjB^ + 



2. for any solution x of Ax — b, B2B2X is k-sparse, i.e., B2B2X e Sfe, 

where T,k = {z : ||z||o < k}; 

3. the 2k-restricted isometry constant 52k for Bi is sufficiently small (for 
example 52k < \/2 — I). 

Then 

{x-xfAix-x)<Ckal{x)i, (27) 

where (Tk{x)i — min^gE^ ||2; — a;||i, x is obtained via the solution of the mini- 
mization problems 

y = argmiiij^llb - Bjy\\i and x = argminj.{||a;||i : Bix = y} (28) 

and the constant Ck only depends on k, the mutual incoherence A4{Bi) and 
\\B2h. 

Proof. By Assumption 1, Ax = b implies that b = BjBix + BJB2X. By 
Assumption 2, it follows that e = BJB2X is fc-sparse. Then by Theorem 1.3 
of [12] we obtain exact recovery, i.e., if Bix = y, then 

y = argmiUyWh- BjyWi. 

The remainder of the proof is then based on Theorems 2.2, 2.4, and Lemma 2.7. 
Since 

A = B'^B = BJBi + BJB2, 

it follows that 

(a; - xf^A{x - x) = \\Bi{x - x)\\l + \\B2ix - x)\\l, (29) 
where x = argmin^{||a;||i : Bix = y}. By Theorem 2.4, we also have 

II -M|2 / II D l|2|| -||2 / 11-^2111^1 fc 2/ x 

11^2(0; - x)\\i < II-B2II2IF - ^Wi < ^ — - f^ifc Wi, 

where Ci.fc only depend on 5k. 

W.l.o.g. we can assume that Bi has unit norm columns. Otherwise instead 
of the linear equation Bix = y, we can consider the following linear equation: 

{B^S)S-^x = y, 

where S = diag( ij^^-p) and Cj is the i-th column of the identity matrix. Then 
the matrix BiS has unit norm column and therefore, by Lemma 2.7 we have 
that 

||Bi(a; - x)\\l < (1 - M{Bi)) \\x - x\\l + M{Bi) \\x - x\\l. 
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By Theorem 2.4, we have ||x — < — |^ a1{x)i and by Theorem 2.2 we have 
— < J. iJ^(a;)i, where Ci^k and C2,fc only depend on Sk- Combining 
these inequaUties with (29), we get 

{x - xfA{x -x)< ((1 - MiB^))^ + M{B,) + 

This conchidos the proof. □ 

Applying Theorem 4.1 to the matrix A = (^B^~"'^)'^ B^~"''^ , where matrix 
QN-n,N _ [^21^22] as in (14), we obtain the corresponding estimate for the 
stiffness matrix. 

Corollary 4.2. Let G M^'^ be a symmetric positive semidefinite matrix 
of rank N - n and A'^x^ = . Let = (^^-".JV^T^Ar-n.w ^ ^^^^ 
factorization of A^ , where B^~'^'^ = [A21 A22]. If the 2k restricted isometry 
constant 52k for B^~'^'^ is sufficiently small (e.g. if 52k < \/2 — then 

(a;^ - xfA'^ix'' - x) < Cfca^(a;^)i, 

where x is obtained via the solution of the minimization problem 

y = argminJ|fe^-B^-"'^^y||i 

and 

X = argmin^{||z||i : B'^^"'^ z = y}, 
and Ck only depends on k and A1(_B^~"'^). 

Proof. Taking Bi = B^""'^, B2 = 0, the proof follows from Theorem 4.1. □ 

Remark 4.3. Equation (27) gives an estimate on the solution of the i'l-mini- 
mization problem. If we assume that ak{x^)i < CNO'k{x°°)i, which means 
that the best approximation of x^ is as good as the best approximation of x°° , 
then Theorem 4.1 shows that the solution that we get from ^i-minimization is 
as good as the best /c-term approximation of x°" , where x°° is the solution of 
original equation A°°a;°° = b°°. 

Remark 4.4. Equation (27) only gives a good bound, if we have 

{Clk - ^)M{B,) < ^(mLx - 1), 

where yUmax is the largest singular value of Bi. Otherwise we may use the direct 
estimate \\Bi{x'^ — x)\\2 < \\Bi\\2\\x'^ — and then apply Theorem 2.2. 

5 Numerical Experiments 

In this section we present some numerical examples. 
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Figure 3: Left: four levels of hat functions on the interval [—1,1]. Right: nonzeros in the 
corresponding stiffness matrix with 11 levels. 



5.1 Example: ID-Poisson Equation 

Let us first demonstrate that £i-minimization can successfully obtain a sparse 
solution. We consider the Poisson equation 

-u" ^ f on n={-l,l), 

with botmdary conditions u(— 1) = = u{l) and 



+ a(x+i)2 l + {a-xy l + a(,T-i)2 l + a2(x-i)2/' 

where a := 100 • tt. The exact solution of this problem is 

u{x) — arctan(Q:(a; + i)) + arctan(— a(a:; + i)) + arctan(Q! • x) 

+ arctan(— a • x) + arctan(a(a; — |)) + arctan(— Q;(a; — j)) 
+ arctan(Q!(a; — i)) + arctan(— Q;(a; — ^)). 

We apply a finite element method [5] and use classical shape functions 



4>{x) 



l-\x\ if - 1 < a; < 1 
otherwise. 



(30) 



The different refinement levels are given by 

(j)kAx) ■■= 2-^/"^ (P{2''-^{x + !)-£), keN, £=!,..., 2'' -1. 

Here, the scaling factor 2^'^'^^ is used to make the diagonal elements of the 
stiffness matrix equal to 1. We then have 



1 



■ifc-i 



{x + l)-i\ if - 1 + 4^ < X < -1 + ^ 



otherwise, 
see the left part of Figure 3 for an illustration. On the A^th level, we then have 



N 



J2i2'' - 1) = 2^+1 ~{N + 2) 



k=l 
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Figure 4: Left: exact solution and approximate solution of the example in Section 5.1 at 
level 8. Right: hat functions that are selected by the £i-minimization problem. 



generating functions and we obtain a stiffness matrix that has a sparsity as 
depicted on the right part of Figure 3. 

For this problem we have numerically computed the mutual incoherence 
and the restricted isometry constant. The numerical results indicate that for 

independently of the size of the matrices, and that 6i — 0, 62 — 0.8165, and for 
all fc > 1, we have (5^ > 1. 

On level 8 we used the matrix of size 255 x 502. The least squares 

solution of A^x = b has 495 nonzeros, while the minimum ^i-solution only 
had 57 nonzeros. The left part of Figure 4 depicts the exact solution and the 
approximate solution at level 8. There is no obvious difference and the relative 
error in ^2-iiorm is 0.0627. The right part of Figure 4 shows that our method 
refines properly at points with large gradients. 



5.2 Application of Algorithm 1 to a ID-Poisson Equation 

To illustrate the behavior of Algorithm 1, wc consider the Poisson equation 

/-„- „"-0 (WnTrfx r, (1007r)^(x-0.5) rPO-f-in 

K-i),z.(i)) = (o,o). 

The exact solution of this problem is 

u — arctan(100 n ■ x) + arctan(— 100 tt) • x 

+ arctan(100 7r(2; — i)) + arctan(— 100 7r)(a; — i). 

We applied four refinement steps of Algorithm 1 starting from level 4. In turns 
out that starting from level 3, a straightforward refinement process does not 
work, since some of the singularities are lost; this is a point where our algorithm 
would need to backtrack, see Section 3.2. 

For the starting level 4, the size of matrix An in (13) is 2'' — (r + 1) = 11 
for r = 4. The size of A21 is (2'' - 1) x (2*^ - (r + 1)) = 15 x 11. The size of A22 
is 15 X 15. 
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Figure 5: Exact and approximate solution for Example 5.2 using Algorithm 1. 



Table 1: Results of Algorithm 1 for Example 5.2. The first column is the refinement step, 
column 2 gives the size of the matrix used for the £i-minimization problem, column 3 gives 
the £o-norm of the £i-minimal solution, column 4 gives the size of the matrix used for the 
classical FEM, column 5 gives the FEM solution of the problem at different levels, and column 
6 gives the ratio of i^o-norms of the £i -solution and the FEM solution. 



step 


size of ^i-niatrix 


ll^llo 


size of FEM matrix 


ll-^llo 


lkllo/||a;||o 


1 


15 X 26 


13 


15 X 15 


15 


0.867 


2 


29 X 42 


23 


31 X 31 


31 


0.742 


3 


59 X 72 


41 


63 X 63 


63 


0.651 


4 


113 X 126 


67 


127 X 127 


127 


0.528 


5 


191 X 204 


103 


255 X 255 


255 


0.404 



At each step we determine the new support using ^i-minimization and then 
refine these nodes and aU necessary higher level nodes according to Algorithm 1 , 
see Figure 5. In Table 1 we present the results of four refinement steps of 
Algorithm 1. 

5.3 Application of Algorithm 1 to a 2D-Poisson Equation 

As a second example for Algorithm 1, we consider the Poisson equation 

-£-0 = ^("'^) [0,l]x[0,l] 

with u{x, y) — Q on the boundary of [0, 1] x [0, 1] where /(x, y) = —20x{x — 1) — 
20y{y — 1). The original solution is u{x,y) — 10xy(x — — 1). 

We use piecewise linear generating functions of the form Ui{x,y) = atx + 
biy + Ci on each triangle. Figure 6 depicts the refinement step from the first 
level (left) to the second level (right). The basis function on level 1 is plotted 
in Figure 7. 

Based on the triangulation on level j and the triangles {A^ }^^-^, we obtain 
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Figure 6: Triangulation on first and second level 






Figure 7: Basis function in the first level 



the generating functions: 



0j,(a,6)(a;,y) = 







a) + 


1 


{x, y) e Aj,fe, 






a)- 


2J(y-6) + l 


{x,y) e Aj,k2 


1 




-b) 


+ 1 


{x,y) e Aj,fe3 


r 


-2i{x 


-a) 


+ 1 


{x,y) e Aj,k^ 




-2i{x 


-a) 


+ 2i{y-b) + l 


(x,y)e Aj^k, 






b) + 


1 


{x, y) e Aj,fcg 



where (a, 6) is the center of the basis function. In general on level n we have 
(2" — 1)^ basis functions. 

We applied four refinement steps of Algorithm 1 starting as first step from 
level 2. At each step we determine the new support using ^i-minimization 
and then refine these nodes and all necessary higher level nodes according to 
Algorithm 1 (see Figure 8). In Table 2 we present the results of four refinement 
steps of Algorithm 1. For the starting level 2, the size of matrix An in (13) is 
Eiri(2'-1)^ Iforr = 2. The size of ^21 is (2''- 1)^ x ^[^^^(2' - 1)^) = 9 x 1. 
The size of A22 is 9 x 9. 

6 Efficiency Estimation 

The algorithm as described in Section 3.2 relies on the solution of a linear 
program (LP) in each iteration. Thus, for it to be successful in practice the 
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Figure 8: Approximate solutions obtained by Algorithm 1 on levels 2,3,4,5,6 and exact 
solution (last) for Example 5.3. 



Table 2: Results of Algorithm 1 for Example 5.3. The first column is the refinement step, 
column 2 gives the size of the matrix used for the £i-minimization problem, column 3 gives 
the ^o-norm of the ^i-minimal solution, column 4 gives the size of the matrix used for the 
classical FEM, column 5 gives the FEM solution of the problem at different level, and column 
6 gives the ratio of £o-iior™s of the ^i-solution and the FEM solution. 



step 


size of i'l-matrix 


ll^llo 


size of FEM matrix 


Mo 


lkllo/||:r||o 


1 


9 X 10 


8 


9x9 


9 


0.889 


2 


43 X 51 


42 


49 X 49 


49 


0.857 


3 


237 X 244 


114 


225 X 225 


225 


0.507 


4 


816 X 823 


172 


961 X 961 


961 


0.179 


5 


1352 X 1359 


190 


3969 X 3969 


3969 


0.048 
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savings in the size of the matrices have to be large enough to compensate for the 
higher solution times of LPs compared to the standard finite element methods. 
Clearly, this depends on a number of factors that are hard to predict: the 
PDE, the data, the required accuracy, the concrete implementation, the usage 
of geometric refinement processes etc. This results in different sizes of the 
refinement tree and different matrices with varying degrees of sparsity. We will 
nevertheless try to get a rough idea of the efficiency. Since our method can be 
seen as a generic way of controlling the refinement process, we compare it against 
the classical finite element method without refinement. We use the examples of 
Section 5 guideline. 

Note that we restrict our attention to the case of solving LPs, although there 
are different methods for compressed sensing that yield similar results as the 
LP-based methods; for instance, orthogonal matching pursuit might be applied, 
see Tropp [41], Donoho, Elad, and Temlyakov [25], Tropp and Gilbert [43] and 
Needell and Tropp [35]. It remains to be seen whether these methods are com- 
petitive for the application discussed in this paper, especially with respect to 
sparse matrices. 

Let us first consider worst case computing times. In the example in Sec- 
tion 5.1, the number of basis functions at level fc is 2*^ — 1, which is also the 
size of the matrix A'' in the equation system (12), which has to be solved by 
a "classical" finite element method. If we use a dense solver this takes 0(2^'^) 
time for each solution. In comparison, dense interior point algorithms for linear 
programming require about 0{n^'^L) time for an LP of dimension n, where L 
is the encoding size of the LP. If the LP is dense, the encoding length includes 
at least one bit for each entry of the matrix and thus is of size at least nm, 
where m is the number of constraints. In our case, we have n = 2^°+^ — {k + 2) 
and m = 2*^ — 1, if we would start our method at level k. Hence, a very optimistic 
estimate of the running time in the dense case would be 0{n^'^m,n) = 0{2^'^^). 

Let us investigate the selection process of the compressed sensing approach 
for this example, see Section 3.2. Assume that we are at iteration k and we have 
the set C^~^ of basis functions with ik-i = #C'^^^. The refined set of basis 
functions C'^ has size at most 3£k-i, because each basis function is subdivided 
into three new basis functions (some of the new basis functions might coincide). 
Assuming that we select at least a fraction of a G (0, 1] among the basis fimctions 
of C'^, the new iteration has at most £k = ^fe-i + a3£k-i = (1 + 3a)£fc_i basis 



Table 3: Running times for solving the LPs of the example in Section 5.1; m and n denote 
the number of rows and columns of the constraint matrix, "time" refers to the running time 
in seconds, and ||z||o gives the number of nonzeros in the solution. 



level 


m 


n 


time 


\\4o 


7 


127 


247 


0.01 


7 


8 


255 


502 


0.03 


8 


9 


511 


1013 


0.10 


9 


10 


1023 


2036 


0.46 


10 


11 


2047 


4083 


1.68 


11 


12 


4095 


8178 


5.74 


12 


13 


8191 


16369 


22.24 


13 


14 


16383 


32752 


88.32 


14 
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120 



80 



40 




I 7 , , , , , 

5000 10000 15000 20000 25000 30000 35000 

Figure 9: Solution times for solving the linear programs of Table 3 versus the number of 
columns of the matrix. 



functions; compare Algorithm 1. Therefore, at iteration k, we have at most 
(1 + 3a)'' = 2'='°S2(i+3a) basis functions. For the compressed sensing approach 
to be successful with respect to the finite element (dense) worst case times above, 
we would need 

log2(l + 3a) < 3.5/5.5 w 0.64 
^ a < i (20 <54 - 1) « 0.19. 

Hence, if we assume that the compressed sensing method yields a reduction 
rate a below approximately 0.19, it should be effective, if we assume worst case 
running times. In the results of Table 1, we have a around 0.5. 

The above estimation is based on the dense worst case running times. Since 
the matrix is sparse, wc can rather assume that the running times for solving the 
equation system are about linear with respect to the size [1, 4, 40]. To estimate 
the running time for linear programming, we have to resort to experiments 
based on the data of Section 5.1. We use the matrices as they would result from 
starting our algorithm at levels 7 to 14. The results are shown in Table 3 and 
Figure 9. We used the barrier solver of CPLEX with additional crossover to 
recover a basic solution. The running times arc with respect to an Intel Core 2 
Quad Core with 2.66 GHz. In Table 3, a seems to be of order a ^ k. Moreover, 
the results suggest a growth of the running time that is lower than quadratic. 
This seems to be a positive sign, which at least does not rule out a possible 
practical effectiveness of our approach. It, however, would require a much more 
thorough computational study to reach definite conclusions. 

Conclusion 

As mentioned in the introduction, many issues of the approach presented in 
this paper have not yet been resolved and many variations arc possible. For 
instance, it is obvious that a similar approach could be derived using other 
dictionaries, e.g., wavelets, instead of finite element functions. Furthermore, 
for practical instances, the solution of the £i-minimization problem becomes 
an issue. One approach would be to apply different algorithms, for instance. 
Orthogonal Matching Pursuit, see [19, 37, 36]. Moreover, the special structure 
of the stiffness matrices can be exploited and techniques adapted to the iterative 
procedure could be developed. 
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